{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from sklearn.datasets import load_iris\n",
    "from sklearn.metrics import accuracy_score\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.preprocessing import StandardScaler"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "class LogisticRegression:\n",
    "    \"\"\" Классификатор на основе логистической регрессии для двух категорий.\n",
    "    \n",
    "    Параметры\n",
    "    ---------\n",
    "    random_state : int\n",
    "      Начальное состояние генератора случаных чисел для первичной инициализации весов\n",
    "    \"\"\"\n",
    "    \n",
    "    def __init__(self, random_state=None):\n",
    "        self._random_state = random_state\n",
    "\n",
    "        if self._random_state:\n",
    "            self.__rand = np.random.RandomState(self._random_state)\n",
    "        else:\n",
    "            self.__rand = np.random\n",
    "    \n",
    "    \n",
    "    def fit(self, X, y, eta=0.01, max_iter=1000, method='GD', cbatch=20, tol=1e-7, sigma=0.1):\n",
    "        \"\"\" Подгонка данных\n",
    "        \n",
    "        Параметры\n",
    "        ---------\n",
    "        X : {array-like}\n",
    "          Матрица обучающих данных\n",
    "        y : array-like\n",
    "          Список категорий\n",
    "        tol : number\n",
    "          Минимальная разница между текущей ошибкой и предыдущей для прерывания оптимизации. По умолчанию 1e-7\n",
    "        sigma : number\n",
    "          Коэфициент импульса. Используется в методе SGD-MOMENTUM, SGD-NESTEROV и SGD-RMSPROP\n",
    "        eta : float\n",
    "          Скорость обучения модели. По умолчанию 0.01\n",
    "        method : str\n",
    "          Метод оптимизации. По умолчанию GD.\n",
    "          Может принимать значения:\n",
    "            GD - градиентный спуск\n",
    "            SGD - стохастический градиентный спуск\n",
    "            SGD-MOMENTUM - стохастический градиентный спуск по стандартной формуле моментов\n",
    "            SGD-NESTEROV - стохастический градиентный спуск по формуле моментов Нестерова\n",
    "            SGD-RMSPROP - стохастический градиентный спуск по формуле RMSPROP\n",
    "        cbatch : int\n",
    "          Количество элементов в каждом батче градиентного спуска. Используется только для методов SGD. По умолчанию 20.\n",
    "        max_iter: int\n",
    "          Количество эпох оптимизации. По умолчанию 1000.\n",
    "        \"\"\"\n",
    "        \n",
    "        method = method.upper()\n",
    "        \n",
    "        method_list = ('GD', 'SGD', 'SGD-MOMENTUM', 'SGD-NESTEROV', 'SGD-RMSPROP')\n",
    "        assert method in method_list, f\"method должен принимать значения {', '.join(method_list)}\"\n",
    "        assert 0<=tol<=1, \"Значение tol должно быть от 0 до 1\"\n",
    "        assert 0<=sigma<=1, \"Значение sigma должно быть от 0 до 1\"\n",
    "        assert 0<=cbatch, \"Значение cbatch должно быть больше или равно 0\"\n",
    "        assert 0<=max_iter, \"Значение max_iter должно быть больше или равно 0\"\n",
    "        assert 0<=eta, \"Значение eta должно быть больше или равно 0\"\n",
    "        assert len(set(y))==2, \"Количество категорий в y должно быть две\"\n",
    "            \n",
    "        # Добавляем x0=1 в каждый обучающий образец\n",
    "        X = np.c_[np.ones(X.shape[0]), X]\n",
    "        \n",
    "        # Инициализируем веса\n",
    "        self._init_weights(X)\n",
    "        \n",
    "        # Переменная с ошибками каждой эпохи\n",
    "        self.cost_ = []\n",
    "        \n",
    "        cost_prev = np.inf\n",
    "        \n",
    "        # Счетчик количества эпох обучения\n",
    "        self.epoch_count_ = 0\n",
    "        \n",
    "        # Оптимизация методом градиентного спуска\n",
    "        if method == 'GD':\n",
    "            for _ in range(int(max_iter)):\n",
    "                \n",
    "                self.epoch_count_ += 1\n",
    "                \n",
    "                cost = self._update_weights_gd(X, y, eta=eta)\n",
    "                \n",
    "                self.cost_.append(cost)\n",
    "                \n",
    "                # Если абсолютное значение разницы текущей ошибки и предыдущей меньше tol, то завершаем обучение\n",
    "                if np.abs(cost_prev - cost) < tol:\n",
    "                    return\n",
    "                \n",
    "                # Текущая ошибка становится предыдущей перед переходом к следующей эпохе\n",
    "                cost_prev = cost\n",
    "                \n",
    "        elif method in ('SGD', 'SGD-MOMENTUM', 'SGD-NESTEROV', 'SGD-RMSPROP'):\n",
    "            \n",
    "            # Инициализируем историю для алгоритмов с моментами и RMSPROP\n",
    "            v_prev = np.zeros(X.shape[1])\n",
    "            \n",
    "            for _ in range(int(max_iter)):\n",
    "                \n",
    "                self.epoch_count_ += 1\n",
    "                \n",
    "                # Перемешиваем обучающие данные\n",
    "                X, y = self._shuffle(X, y)\n",
    "                \n",
    "                # Получаем разделенные значение X и y с количеством елементов в блоке равным cbatch\n",
    "                X_batches, y_batches = self._make_splitted_batches(X=X, y=y, cbatch=cbatch)\n",
    "                \n",
    "                # Переменная для накопления ошибок при расчете каждого батча\n",
    "                epoch_cost = 0 \n",
    "                \n",
    "                # Для каждого батча производим оптимизацию и расчет ошибки в зависимости от типа SGD\n",
    "                for X_batch, y_batch in zip(X_batches, y_batches):\n",
    "                    if method == 'SGD':\n",
    "                        epoch_cost += self._update_weights_gd(X_batch, y_batch, eta=eta)\n",
    "                    elif method == 'SGD-MOMENTUM':\n",
    "                        cost, v_prev = self._update_weights_momentum(X_batch, y_batch, eta=eta, v_prev=v_prev, sigma=sigma)\n",
    "                        epoch_cost += cost\n",
    "                    elif method == 'SGD-NESTEROV':\n",
    "                        cost, v_prev = self._update_weights_nesterov(X_batch, y_batch, eta=eta, v_prev=v_prev, sigma=sigma)\n",
    "                        epoch_cost += cost\n",
    "                    elif method == 'SGD-RMSPROP':\n",
    "                        cost, v_prev = self._update_weights_rmsprop(X_batch, y_batch, eta=eta, v_prev=v_prev, sigma=sigma)\n",
    "                        epoch_cost += cost\n",
    "                \n",
    "                # Считаем среднюю ошибку в батчах для эпохи\n",
    "                cost = epoch_cost/len(X_batches)\n",
    "                self.cost_.append(cost)\n",
    "                \n",
    "                # Если абсолютное значение разницы текущей ошибки и предыдущей меньше tol, то завершаем обучение\n",
    "                if np.abs(cost_prev - cost) < tol:\n",
    "                    return\n",
    "                \n",
    "                # Текущая ошибка становится предыдущей перед переходом к следующей эпохе\n",
    "                cost_prev = cost \n",
    "    \n",
    "    \n",
    "    def _update_weights_gd(self, X, y, eta):\n",
    "        \"\"\" Обновление весов методом градиентного спуска\n",
    "        \n",
    "        Параметры\n",
    "        ---------\n",
    "        X : {array-like}\n",
    "          Матрица обучающих данных\n",
    "        y : array-like\n",
    "          Список категорий\n",
    "        eta : float\n",
    "          Скорость обучения модели.\n",
    "          \n",
    "        Результат\n",
    "        ---------\n",
    "        cost : float\n",
    "          Размер ошибки для логистической регрессии\n",
    "        \"\"\"\n",
    "        \n",
    "        # Вычисляем функцию общего входа для каждого обучающего элемента\n",
    "        net_input = self._net_input(X=X, w=self._w)\n",
    "        \n",
    "        # Вычисляем значение функции активации (сигмоида) для каждого обучающего элемента\n",
    "        output = self._activation(net_input)\n",
    "        \n",
    "        # Обновляем веса по формуле градиента логистичесикй регрессии\n",
    "        self._w += eta * X.T.dot(y-output)\n",
    "        \n",
    "        # Рассчитываем текущую ошибку\n",
    "        cost = self._cost_calc(z=output, y=y)\n",
    "        \n",
    "        return cost  \n",
    "    \n",
    "    \n",
    "    def _update_weights_rmsprop(self, X, y, eta, v_prev, sigma=0.9):\n",
    "        \"\"\" Обновление весов методом RMSPROP\n",
    "        \n",
    "        Параметры\n",
    "        ---------\n",
    "        X : {array-like}\n",
    "          Матрица обучающих данных\n",
    "        y : array-like\n",
    "          Список категорий\n",
    "        eta : float\n",
    "          Скорость обучения модели.\n",
    "        sigma : number\n",
    "          Коэфициент импульса формулы оптизизации rmsprop\n",
    "          sigma * v_prev + (1 - sigma) * np.dot(gradient, gradient)\n",
    "        v_prev : array-like\n",
    "          Предыдущее значение формулы sigma * v_prev + (1 - sigma) * np.dot(gradient, gradient)\n",
    "          \n",
    "        Результат\n",
    "        ---------\n",
    "        cost, v_current : tuple\n",
    "          Размер ошибки для логистической регрессии \n",
    "          и текущее значение формулы sigma * v_prev + (1 - sigma) * np.dot(gradient, gradient)\n",
    "        \"\"\"\n",
    "        \n",
    "        self._X = X\n",
    "        eps = np.finfo(np.float64).eps\n",
    "        \n",
    "        # Вычисляем функцию общего входа для каждого обучающего элемента\n",
    "        net_input = self._net_input(X=X, w=self._w)\n",
    "        # Вычисляем значение функции активации (сигмоида) для каждого обучающего элемента\n",
    "        output = self._activation(net_input)\n",
    "        \n",
    "        # Вычисляем градиент\n",
    "        gradient = X.T.dot(y-output)\n",
    "        self._gradient = gradient\n",
    "        \n",
    "        v_current = sigma * v_prev + (1 - sigma) * np.dot(gradient, gradient)\n",
    "        # Обновляем веса по формуле RMSPROP\n",
    "        self._w += (eta/np.sqrt(v_current + eps))*gradient\n",
    "        \n",
    "        # Рассчитываем текущую ошибку\n",
    "        cost = self._cost_calc(z=output, y=y)\n",
    "        \n",
    "        return cost, v_current\n",
    "    \n",
    "    \n",
    "    def _update_weights_momentum(self, X, y, eta, v_prev, sigma=0.9):\n",
    "        \"\"\" Обновление весов методом градиентного спуска\n",
    "        \n",
    "        Параметры\n",
    "        ---------\n",
    "        X : {array-like}\n",
    "          Матрица обучающих данных\n",
    "        y : array-like\n",
    "          Список категорий\n",
    "        eta : float\n",
    "          Скорость обучения модели.\n",
    "        sigma : number\n",
    "          Коэфициент импульса формулы оптизизации моментов\n",
    "          sigma * v_prev + eta * X.T.dot(y-output)\n",
    "        v_prev : array-like\n",
    "          Предыдущее значение формулы sigma * v_prev + eta * X.T.dot(y-output)\n",
    "          \n",
    "        Результат\n",
    "        ---------\n",
    "        cost, v_current : tuple\n",
    "          Размер ошибки для логистической регрессии и текущее значение формулы sigma * v_prev + eta * X.T.dot(y-output)\n",
    "        \"\"\"\n",
    "        \n",
    "        # Вычисляем функцию общего входа для каждого обучающего элемента\n",
    "        net_input = self._net_input(X=X, w=self._w)\n",
    "        # Вычисляем значение функции активации (сигмоида) для каждого обучающего элемента\n",
    "        output = self._activation(net_input)\n",
    "        \n",
    "        v_current = sigma * v_prev + eta * X.T.dot(y-output)\n",
    "        # Обновляем веса по формуле моментов\n",
    "        self._w += v_current\n",
    "        \n",
    "        # Рассчитываем текущую ошибку\n",
    "        cost = self._cost_calc(z=output, y=y)\n",
    "        \n",
    "        return cost, v_current\n",
    "    \n",
    "    \n",
    "    def _update_weights_nesterov(self, X, y, eta, v_prev, sigma=0.9):\n",
    "        \"\"\" Обновление весов методом градиентного спуска\n",
    "        \n",
    "        Параметры\n",
    "        ---------\n",
    "        X : {array-like}\n",
    "          Матрица обучающих данных\n",
    "        y : array-like\n",
    "          Список категорий\n",
    "        eta : float\n",
    "          Скорость обучения модели.\n",
    "        v_prev :  array-like\n",
    "          Предыдущее значение формулы sigma*v_prev - eta * X.T.dot(y-output)\n",
    "        sigma : number\n",
    "          Коэффициент импульса формулы моментов Нестерова\n",
    "          \n",
    "        Результат\n",
    "        ---------\n",
    "        cost, v_current : tuple\n",
    "          Размер ошибки для логистической регрессии и текущее вычисление sigma*v_prev - eta * X.T.dot(y-output)\n",
    "        \"\"\"\n",
    "        \n",
    "        # Вычисляем веса для места ближе к точке, где окажемся на след.шаге\n",
    "        temp_w = self._w - sigma*v_prev\n",
    "        \n",
    "        # Вычисляем функцию общего входа для каждого обучающего элемента с модифицированными весами\n",
    "        net_input = self._net_input(X=X, w=temp_w)\n",
    "        # Вычисляем значение функции активации (сигмоида) для каждого обучающего элемента с модифицированными весами\n",
    "        output = self._activation(net_input)\n",
    "        \n",
    "        v_current = sigma*v_prev - eta * X.T.dot(y-output)\n",
    "        \n",
    "        # Обновляем веса по формуле Нестерова\n",
    "        self._w -= v_current\n",
    "        \n",
    "        # Для каждой эпохи расчитываем ошибку\n",
    "        cost = self._cost_calc(z=output, y=y)\n",
    "        \n",
    "        return cost, v_current\n",
    "    \n",
    "    \n",
    "    def _make_splitted_batches(self, X, y, cbatch):\n",
    "        \"\"\" Разделение обучающих данных на батчи для градиентного спуска\n",
    "        \n",
    "        Параметры\n",
    "        ---------\n",
    "        X : {array-like}\n",
    "          Матрица обучающих данных\n",
    "        y : array-like\n",
    "          Список категорий\n",
    "        cbatch : int\n",
    "          Количество элементов в каждом батче\n",
    "          \n",
    "        Результат\n",
    "        ---------\n",
    "        X_batches, y_batches: (list, list)\n",
    "          Списки батчей для X и y\n",
    "        \"\"\"\n",
    "        \n",
    "        length_X = len(X)\n",
    "        assert cbatch <= length_X, 'cbatch не должен превышать количества строк в X'\n",
    "\n",
    "        # Блок разделения данных обучения на части. В каждом батче количество обучающих данных равно cbatch\n",
    "        start_index = 0\n",
    "        end_index = cbatch\n",
    "        X_batches =[]\n",
    "        y_batches = []\n",
    "        for _ in range(length_X//cbatch):\n",
    "            X_batches.append(X[start_index:end_index])\n",
    "            y_batches.append(y[start_index:end_index])\n",
    "            start_index += cbatch\n",
    "            end_index += cbatch\n",
    "        # Добавляем оставшиеся элементы, которые не были добавлены в цикле\n",
    "        if len(X[start_index:]) > 0:\n",
    "            X_batches.append(X[start_index:])\n",
    "            y_batches.append(y[start_index:])\n",
    "        \n",
    "        return X_batches, y_batches\n",
    "        \n",
    "\n",
    "    def _cost_calc(self, z, y):\n",
    "        \"\"\" Расчет ошибки для логистической регрессии\n",
    "        \n",
    "        Параметры\n",
    "        ---------\n",
    "        z : array-like or number\n",
    "          Значение функции общего входа. Список или число.\n",
    "        y : array-like\n",
    "          Значение корректных категорий. Список или число.\n",
    "          \n",
    "        Результат\n",
    "        ---------\n",
    "        cost : number\n",
    "          Ошибка логистической регрессии\n",
    "        \"\"\"\n",
    "  \n",
    "        return -np.dot(y, np.log(z)) - np.dot(1-y, np.log(1-z))\n",
    "    \n",
    "    \n",
    "    def _init_weights(self, X):\n",
    "        \"\"\" Инициализация весов\n",
    "        \"\"\"\n",
    "        # Инициализируем веса близкими к 0 значениями\n",
    "        self._w = self.__rand.normal(loc=0, scale=0.01, size=X.shape[1])\n",
    "      \n",
    "    \n",
    "    def _shuffle(self, X, y):\n",
    "        \"\"\" Перемешивание обучающих данных\n",
    "        \n",
    "        Параметры\n",
    "        ---------\n",
    "        X : {array-like}\n",
    "          Матрица обучающих данных\n",
    "        y : array-like\n",
    "          Список категорий\n",
    "          \n",
    "        Результат\n",
    "        ---------\n",
    "        X, y: tuple\n",
    "          Перемешанные X и соответствующие им значения y\n",
    "        \"\"\"\n",
    "        r = self.__rand.permutation(len(X))\n",
    "        return X[r], y[r]\n",
    "      \n",
    "        \n",
    "    def _net_input(self, X, w=None):\n",
    "        \"\"\" Вычисляет общий вход z\n",
    "        \n",
    "        Параметры\n",
    "        ---------\n",
    "        X : {array-like}\n",
    "          Массив признаков\n",
    "        w : array-like\n",
    "          Веса. Если w = None, то будут взяты веса из self._w\n",
    "          \n",
    "        Результат\n",
    "        ---------\n",
    "        z : float\n",
    "          Значение общего входа\n",
    "        \"\"\"\n",
    "        if w is None:\n",
    "            w = self._w\n",
    "            \n",
    "        return np.dot(X, w)\n",
    "    \n",
    "    \n",
    "    def _activation(self, z):\n",
    "        \"\"\" Выполняет сигмоидальную активацию\n",
    "        \n",
    "        Параметры\n",
    "        ---------\n",
    "        z : float\n",
    "          значение общего входа\n",
    "          \n",
    "        Результат\n",
    "        ---------\n",
    "        f : float\n",
    "          Значение сигмоидальной функции активации\n",
    "        \"\"\"\n",
    "        # Ограничиваю сверху размером 36, так как при большем значении функция дает результат 1\n",
    "        # и из-за этого невозможно расчитать ошибку, так как выражение np.log(1-1) выдает исключение\n",
    "        return 1 / (1 + np.exp(-np.clip(z, -250, 36)))\n",
    "    \n",
    "    \n",
    "    def predict(self, X):\n",
    "        \"\"\" Предсказание категории\n",
    "        \n",
    "        Параметры\n",
    "        ---------\n",
    "        X : {array-like}\n",
    "          Массив признаков\n",
    "          \n",
    "        Результат\n",
    "        ---------\n",
    "        predict : array-like\n",
    "          Массив категорий для массива признаков\n",
    "        \"\"\"\n",
    "        # Добавляем x0=1 в каждый обучающий образец\n",
    "        X = np.c_[np.ones(X.shape[0]), X]\n",
    "        \n",
    "        return np.where(self._net_input(X) > 0.0, 1, 0) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "X, y = load_iris(return_X_y=True)\n",
    "# Оставляем только признаки 0 и 1\n",
    "X = X[:100]\n",
    "y = y[:100]\n",
    "\n",
    "X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3)\n",
    "scaler = StandardScaler()\n",
    "scaler.fit(X_train)\n",
    "X_train_std = scaler.transform(X_train)\n",
    "X_test_std = scaler.transform(X_test)\n",
    "\n",
    "lr = LogisticRegression(random_state=42)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Метод градиентного спуска\n",
      "accuracy: 1.0\n",
      "Количество эпох: 184\n"
     ]
    }
   ],
   "source": [
    "print(\"Метод градиентного спуска\")\n",
    "lr.fit(X_train, y_train, method='GD', tol=1e-3)\n",
    "y_predict = lr.predict(X_test)\n",
    "print('accuracy:', accuracy_score(y_test, y_predict))\n",
    "print('Количество эпох:', lr.epoch_count_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Метод стохастического градиентного спуска\n",
      "accuracy: 1.0\n",
      "Количество эпох: 32\n"
     ]
    }
   ],
   "source": [
    "print(\"Метод стохастического градиентного спуска\")\n",
    "lr.fit(X_train, y_train, method='SGD', cbatch=1, tol=1e-3)\n",
    "y_predict = lr.predict(X_test)\n",
    "print('accuracy:', accuracy_score(y_test, y_predict))\n",
    "print('Количество эпох:', lr.epoch_count_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Метод стохастического градиентного спуска по стандартной формуле моментов\n",
      "accuracy: 1.0\n",
      "Количество эпох: 10\n"
     ]
    }
   ],
   "source": [
    "print(\"Метод стохастического градиентного спуска по стандартной формуле моментов\")\n",
    "lr.fit(X_train, y_train, method='SGD-MOMENTUM', cbatch=1, tol=1e-3, sigma=0.9)\n",
    "y_predict = lr.predict(X_test)\n",
    "print('accuracy:', accuracy_score(y_test, y_predict))\n",
    "print('Количество эпох:', lr.epoch_count_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Метод стохастического градиентного спуска по формуле моментов Нестерова\n",
      "accuracy: 1.0\n",
      "Количество эпох: 11\n"
     ]
    }
   ],
   "source": [
    "print(\"Метод стохастического градиентного спуска по формуле моментов Нестерова\")\n",
    "lr.fit(X_train, y_train, method='SGD-NESTEROV', cbatch=1, tol=1e-3, sigma=0.9)\n",
    "y_predict = lr.predict(X_test)\n",
    "print('accuracy:', accuracy_score(y_test, y_predict))\n",
    "print('Количество эпох:', lr.epoch_count_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Метод стохастического градиентного спуска по формуле RMSPROP\n",
      "accuracy: 1.0\n",
      "Количество эпох: 32\n"
     ]
    }
   ],
   "source": [
    "print(\"Метод стохастического градиентного спуска по формуле RMSPROP\")\n",
    "lr.fit(X_train, y_train, method='SGD-RMSPROP', cbatch=1, tol=1e-3, sigma=0.9)\n",
    "y_predict = lr.predict(X_test)\n",
    "print('accuracy:', accuracy_score(y_test, y_predict))\n",
    "print('Количество эпох:', lr.epoch_count_)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
